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Abstract 

Trapped Bose atoms cooled down to temperatures below the Bose-Einstein con- 
densation temperature are considered. Stationary solutions to the Gross-Pitaevskii 
equation (GPE) define the topological coherent modes, representing nonground- 
state Bose-Einstein condensates. These modes can be generated by means of alter- 
nating fields whose frequencies are in resonance with the transition frequencies be- 
tween two collective energy levels corresponding to two different topological modes. 
The theory of resonant generation of these modes is generalized in several aspects: 
Multiple-mode formation is described; a shape-conservation criterion is derived, im- 
posing restrictions on the admissible spatial dependence of resonant fields; evolution 
equations for the case of three coherent modes are investigated; the complete sta- 
bility analysis is accomplished; the effects of harmonic generation and parametric 
conversion for the topological coherent modes are predicted. All considerations are 
realized both by employing approximate analytical methods as well as by numeri- 
cally solving the GPE. Numerical solutions confirm all conclusions following from 
analytical methods. 



PACS numbers: 03.75.-b, 03.75.Fi, 05.45.-a 
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1 Introduction 



Dilute Bose gases at low temperatures, when almost all atoms are in Bose-Einstein con- 
densate, are well described by the GPE (see reviews [I] [H]). Since the latter is a type of 
the nonlinear Schrodinger equation, it should possess the whole spectrum of stationary 
solutions. In the presence of a trapping potential, the related energy spectrum is, in gen- 
eral, discrete. Stationary solutions to the GPE are, by definition, the topological coherent 
modes, whose ground state describes the standard Bose-Einstein condensate, while the 
higher states correspond to nonground-state Bose-Einstein condensates 0. 

It is worth saying a few words recalling where the name topological coherent modes 
comes from. Different stationary solutions to the GPE, associated with different energy 
levels, display principally different spatial shapes, in particular, different number of zeros. 
Because of their distinct spatial topology, the modes are termed topological. These should 
not be confused with elementary collective excitations, defined by linear Bogolubov - 
De Gennes equations. Such elementary excitations describe small oscillations around a 
given stationary solution and do not change the topology of the latter. Since the GPE is 
nonlinear, its solutions can also be named nonlinear modes, which would stress their prin- 
cipal difference from elementary excitations satisfying the linear Bogolubov - De Gennes 
equations. However, the elementary excitations, produced by strong perturbations, are 
sometimes also called nonlinear. Therefore the term topological, characterizing dissimilar 
modes, seems to be more precise. The topological modes are specified as coherent due to 
the fact that the GPE can be interpreted as an exact equation for coherent states |H] . 

The possibility of resonant generation of arbitrary topological coherent modes was ad- 
vanced in Ref. A particular case of vortex creation was suggested in EH EI]- The 
properties of such modes were also studied in [T3] - [22] and a dipole topological mode was 
excited in experiment [2B|. Bose-Einstein condensates with topological coherent modes 
exhibit a variety of interesting features which could find many applications. Among these 
features, we could mention the following. Mode locking: this is the effect under which 
the fractional mode populations are locked to stay in the vicinity of their initial values 
|HllHll22j- Mathematically, this effect is analogous to self-trapping occurring for atoms 
in double- well potentials [23 123 l2Hj- Critical dynamics: an abrupt qualitative change 
of dynamics of fractional mode populations under an infinitesimally small variation of 
the pumping parameters (HI El HU 122] • A mathematically similar effect in the case of 
double-well potentials is the dynamic phase transition between the Rabi and Josephson 
regimes [2UI23I2ZI- Interference patterns: specific interference fringes arising because of 
differing spatial shapes of different topological modes [HI IS]- Interference current: fastly 
oscillating current, existing even inside a single-well trap [El IS]- Such a type of current, 
arising between two interpenetrating populations, not separated by any barrier, is, on 
occasion, called the internal Josephson effect [2E1 122] • Atomic squeezing: narrowing of 
the dispersion corresponding to the mode population difference, compared to the disper- 
sion related to dipole transitions jHj- A similar feature is illustrated by two-component 
condensates [3U]. Irreducible modes: these are the topological modes that have no lin- 
ear counterparts, so that they cannot be considered as analytical continuations, under 
increasing of nonlinearity, of the related linear modes jTS| EI] • Such modes, in addition 
to strong nonlinearity, require the presence of multiwell potentials [HH 1221 HUH Oil • An 
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investigation of the mode spectrum in the case of single-well traps shows that in this case 
nonlinear modes are reducible and can be treated as analytical continuations of linear 
counterparts IH1 EES E51 ESI • 

In the present paper, we generalize the theory of resonant formation of topological 
coherent modes and study the features that have not been considered in previous publi- 
cations. The most important novel points follows: 

1. The possibility of resonant generation of multiple topological coherent modes is 
described. This can be achieved by subjecting a trapped Bose-Einstein condensate 
to the action of several alternating fields, whose frequencies are tuned to distinct 
transition frequencies related to different modes (Sec. 

2. A general criterion is derived, showing when nonlinear modes cannot be generated 
even if the applied alternating field is in perfect resonance with the corresponding 
transition frequency. This condition relates the spatial dependence of the trapping 
potential and that of the alternating field (Sec. EJ). 

3. Simultaneous generation of two excited coherent modes is studied in detail. The 
resulting dynamical system describes three coexisting nonlinear modes (Sec.EJ). 

4. The phase portrait of the three-mode dynamical system is investigated. All fixed 
points are found and their stability is analysed (Sec. EJ). 

5. Harmonic generation of topological modes is shown to exist. This effect is analogous 
to optical harmonic generation (Sec.[7|). 

6. Parametric conversion for topological modes is another effect having its optical 
counterpart. To realize this effect, it is necessary to subject trapped atoms to the 
action of several alternating fields (Sec. |7j). 

7. The principal feature of this article is that all consideration has been done in two 
ways, by applying approximate analytical methods and also by numerically solving 
the GPE. Numerical solutions confirm all effects described analytically. 



Dilute Bose-condensed gases at low temperature are characterized by a coherent 

field (p(r, t), which is a wave function satisfying the GPE 
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Resonant Generation of Multiple Modes 





where the nonlinear Hamiltonian 
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contains a trapping potential U(r) and the interaction intensity 



m 



(3) 



with m being atomic mass; a s , scattering length; and N, the total number of atoms. The 
wave function is normalized to unity, so that \\<p\\ = 1. The confining potential can be 
modulated by applying an additional field V = V(r, t). 

The topological coherent modes are the solutions to the stationary GPE 



H[(p n ) ip n (r) = E n (p n (r) 



(4) 



where n is a labelling multi-index. The transition frequencies for two distinct modes 
m 7^ n are 

&mn = T" (E m — E n ) . (5) 

n 

The trap modulation is resonant if the frequency of an applied alternating field is tuned 
close to one of the transition frequencies (j5J). This resonant field can induce transitions 
between the considered modes. Everywhere in what follows, talking about resonance, we 
keep in mind the resonant transitions between distinct topological modes. This should 
not be confused with parametric resonance, when one considers a single mode, whose 
width can become divergent under special conditions on the amplitude of a perturbing 
field |SZ|. 

Previously, the resonant excitation of topological modes has been considered for the 
case of a sole resonant field coupling two chosen modes IE] • Now, we turn to the most 
general case, when there are several alternating fields, so that the modulating potential 
is a sum 



K(r,t)=^[V 7 -(r) cosK-t) 



V?(r) 



sm(u)jt) 



(6) 



of fields with different frequencies Uj. Generally, the amplitudes Vj(r) and V-(r) can also 
be slow functions of time, whose temporal variation is slow compared to that of cos^-t), 
such that 



hujj 



Of 



<C 1 



dV< 



Of 



< 1 



(7) 



If these amplitudes are such slow functions of time, then, varying them adiabatically, one 
could induce the Landau-Zener tunneling [2U EH] ■ But here we shall consider only the 
transitions caused by the fastly oscillating cos(ojjt). Each frequency Uj is assumed to be 
close to one of the transition frequencies (0), so that the resonance conditions 



UJ r , 



are valid. 

The modulating potential (0) can be presented as 



V(r,t) = W[Bj(r) e**' + 2$(r) e 



(9) 
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with 

BM^VM-iVjir). (10) 

The modulation can be realized by varying the trapping magnetic fields, by invoking laser 
spoons, or resorting to other means (HI El El HU| ITT]. 

There exist two characteristic quantities, called transition amplitudes; one is the ma- 
trix element 

a mn = N—^ (\ipm\ 2 , 2\ip n \ 2 - \ip m \ 2 ) , (11) 



h 

due to the interatomic interactions (J3J), and another one is the matrix element 



Ann = T {<Pm, B j<fn) , (12) 



h 

related to the amplitude Bj = Bj(r) of the modulating potential ©• Here, (., .) denotes 
the usual scalar product of Schrodinger theory. To avoid intensive power broadening, 
these amplitudes ((UJ and (]T2]) have to satisfy the inequalities 



a. 



OJr. 



< 1 



A 



OJr. 



« 1 , (13) 



where m 7^ n and whose meaning was explained in detail in [Hj. Conditions (|13jl allow for 
an effective generation of nonlinear modes by resonant fields. The first of these restrictions, 
briefly speaking, can be reduced to the limitation on the number of atoms that can be 
transferred to an excited coherent mode jH]. This limiting number of atoms is close to 
the critical number of atoms with attractive interactions, for which the Bose-Einstein 
condensate preserves its stability [HI IH1 EM ESI- Therefore the resonant generation of 
nonlinear modes is feasible for atoms with positive as well as negative scattering lengths. 

The topological coherent modes, being the solutions to the nonlinear eigenproblem (|4*]l . 
do not compulsory form a complete orthonormal basis. However, the modes <£>n(r) can 
always be normalized, so that ||<^ n || — !• And the set {(p n (r)} of all linearly independent 
functions forms a total basis jH], permitting one to look for the solution of the Gross- 
Pitaevskii equation (JTJ) in terms of the presentation 

cp{r,t) = S ^c n (t) y? n (r) exp (- - E n tj , (14) 

where c n (t) are unknown functions of time. Note that for some nonlinear eigenproblems, 
it has been rigorously proved (JOl El 02] that the set of the corresponding eigenfunctions 
forms a complete basis. In our case, it is sufficient to require that the functions c n (t) are 
slowly varying, such that 

1 



OJr. 



dt 



« 1 , (15) 



and conditions (|TH|) hold true. Then the expansion ()14jl is uniquely defined by means of 
the averaging technique j33]. 

It is worth emphasizing that the presentation (14) ideally suits for analysing resonant 
transitions between the coherent topological modes. And it is solely these modes that are 
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the subject of the present paper. We shall not consider here other possible excitations 
that could be produced by nonresonant driving fields and studied by means of the known 
rescaling procedure. Nontopological nonresonant breathing-type oscillations have been 
considered by many authors in the early days of the BEC research (see reviews PQ [Ej)- 
Therefore there is no reason of extending the paper by repeating similar results. This is 
why here we limit ourselves by treating only the resonant generation of topological modes, 
which have not been studied earlier. 

Aiming at exciting particular modes, one should keep in mind, in addition to the 
resonance conditions (JBJ), the symmetry properties of the corresponding wave functions, 
for which the modulating field has to be such that the transition amplitudes (fT2|) be 
nonzero. However, even if all above mentioned conditions hold true, there exists a rather 
general situation when the generation of modes is impossible. 

3 No-Go Theorem for Mode Generation 

It may happen that the applied modulating field is not able to generate higher nonlinear 
modes, but can only lead to the oscillation of the initial wave function, without changing 
its shape. More precisely, let us start with a wave function (p(r,0). After a modulating 
field V(r, t) begins acting on atoms, the function tp(r, 0) is transferred to a function ip(r, t). 
When the shape- conservation condition 

|^(r,t)| = |^(r-a,0)|, (16) 

holds true for time dependent a = a(t), then the shape of the atomic cloud does not change 
in time, but the cloud oscillates as a whole, with its center of mass moving according to 
the dependence a(t). Hence, if we start with a mode (p(r, 0) it can never be transferred 
to another mode. 

By definition, the initial function <^(r,0) = (fo{ r ) presents a nonlinear mode if it 
satisfies the stationary Gross-Pitaevskii equation 

if Mr)] v?o(r) = E p (r) , (17) 

with the nonlinear Hamiltonian (J2J). We assume that this initial mode is real- valued, i.e., 

<p(r,0) = <p (r) = <p*(r) . (18) 

An example for this would be the ground-state of a Bose-Einstein condensate. In addition, 
we are focusing on trapped atoms, which implies that the confining potential U (r) increases 
towards infinity for r = |r| — > oo. Therefore, the trapping condition 

lim (p(r,t) = (19) 

r^oo 

is valid for all t > 0. 

Theorem. Suppose that atoms in a trapping potential U(r), being initially in a real 
mode ifio{r), are subject to the action of a modulating field V(r,t), so that conditions 
(JT7j) to (|19|) are valid. Then the solution of the temporal Gross-Pitaevskii equation (0) 
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preserves the shape of the initial mode, satisfying condition (j!6|) . if and only if the trapping 
potential is harmonic, 

U(r) = A + A 1 -r + Y,A a pr a r?, (20) 

a/3 

where a and (3 are the Cartesian indices, while the modulating field is linear with respect 
to the spatial variables, 

V(r,t) = B (t)+B 1 (t)-T, (21) 

Bo(t) and Bi(i) being arbitrary functions of time. And the center-of-mass motion a = a(t) 
is described by the equation 

m ^- + J2( A ^ + A ^) a + B?(t) = O. (22) 

Proof. Assume that the shape- conservation condition ()16|) holds true, and let us show 
that then the trapping and modulating potentials are necessarily such as in Eqs. (|2U|) and 
(|'2ip. For this purpose, it is useful to invoke the presentation 

<p(r,t) = \ V (r,t)\ exp{zS(r,t)} , (23) 

where S(r, t) is a real action defining the velocity 

v(r,t) = — VS(r,t). (24) 
m 

Substituting the presentation (f23j) into the GPE (JJJ yields PQ-jH] the continuity equation 

|m 2 + V(M 2 v) = (25) 

and the velocity equation 

m = -W e// , (26) 

in which the effective potential 

+ (27) 

Because of the shape-conservation condition (|16p. one has 

9|<p| . da 

dt |y| dt 

Using this in the continuity equation ()25|) results in 

M 2 (v-f 

V dt , 
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from where 

da c 

at \(p\ z 

If here the function c(t) is not zero, then from the trapping condition (JT5J) it follows that 
v — > oo as r — ► oo hence, the action 5 — > oo as r — > oo. But then the limit r — > oo of the 
function (|2*5|) is not defined. Therefore, c(i) = 0. Thus, the velocity (|2*4*|) becomes 

v=|=v( t ), (28) 

which is a function of time only. But then the velocity equation ()26|) tells us that VU e ff 
is also a function of time only, hence the effective potential is linear in r, having the form 

U eff = D (jt)+D 1 {t)-r. (29) 

Taking into account that the initial mode is real- valued, as is conditioned by Eq. (|18|). 
and employing the shape- conservation condition JTBJ), one gets 

|^(r,t)|=^ (r-a) . (30) 

Then the effective potential (|27|) can be written as 

2 

U eff = U(t) - U{r - a) + V{r,t) + E + . (31) 

Since, according to Eq. (J2HJ), the effective potential is linear in r, E is a constant, and 
v = v(t) is a function of time, then the sum of the first three terms in Eq. (j3*Tj) should 
also be linear in r, such that 

U(r) - U(r - a) + V(r, t) = U (t) + Ui(t) • r . (32) 

Thence the effective potential (|31|) takes the form 

2 

U eff = E + + U (t) + Ui(t) ■ r . (33) 

Comparing Eqs. and (|3^j) . we have 

L> (t) = £o + ^ + I7o(*) , Dx(t) = U t (t) • (34) 

Equality can be satisfied only if the trapping and modulating potentials are given by 
Eqs. d23l and d2U). Substituting the latter in Eq. we find 

U (t) = B (t) + a • A - A a p a a a? , 

U?(t) = B?(t) + ( A «P + A Pa) o? . (35) 
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Combining Eqs. (J25j) . and we obtain Eq. (J22"j) for the center-of-mass 

motion. 

Now let us show that Eqs. (|2(J|) and (|21|) are sufficient for the validity of the shape- 
conservation condition (fTHj). Equations an d (ESI)? with the effective potential (|2"7|) 
and with the initial conditions 

|^(r,0)| = ^(r), v(r,0)=0, 

are equivalent to the GPE (JTJ with the initial condition (|18|) . The latter equation is a 
nonlinear Schrodinger equation, which, being complimented by the boundary condition 
JIHJ), possesses a unique solution Hence, Eqs. (J25J) and (ffitj) . with the same boundary 
condition (fTUjl. enjoy a unique solution. These equations, under conditions (p2"Uj) and (JUJ, 
do have a solution satisfying the shape-conservation condition which, according to 
the aforesaid, is a unique solution. This concludes the proof. 

One should not confuse the shape-conservation criterion derived here with the known 
result of the decoupled center-of-mass motion in a harmonic potential. The criterion of 
shape conservation shows when the wave function retains its shape under the action of an 
external alternating field and when the shape is not preserved. The center-of-mass motion 
is a trivial byproduct of our theorem. The principal point is the shape conservation. The 
derived criterion shows that the trapping potential may be harmonic, but, if the driving 
field is nonlinear, then, irrespectively of the center-of-mass motion, the shape of the wave 
function will not be conserved. 

This theorem teaches us that if the trapping potential is harmonic, which is a standard 
situation, then the modulating field, being linear in r, is not able to generate topological 
modes, even if its alternating temporal parts oscillate in an exact resonance with the 
corresponding transition frequencies. It is in full agreement with a nonlinear Ehrenfest- 
Theorem for the mean position and variance of the collective wavefunction 11 . To 
generate such modes, at least one of the conditions (J2Uj) or (J2T|l has to be broken. But 
if a resonant alternating field is linear and the trapping potential is harmonic, then the 
real initial mode just moves in space, without changing its shape. We have checked this 
conclusion by numerically solving the GPE under conditions (|2(J|) and (J21|) and found a 
perfect agreement with the theorem. 

4 Coupling of Three Nonlinear Modes 

Several modes can be generated by applying a modulating field © containing several 
corresponding resonant terms. In general, two cases are admissible: when the excited 
modes are not coupled to each other and when they are coupled. In the former case, 
the overall dynamics consists of the motion of several pairs of modes, each pair being 
separated in its motion from other modes. The motion of such separate pairs of resonant 
modes has been studied in detail earlier [HI El 13 EE HH UH UH EQ 122] • Therefore we now 
concentrate on the case of several coupled modes. 

We consider the case of three coupled modes, whose transition frequencies (j5J) are 
enumerated as u>2i, ct>3i, and U32, keeping in mind that the related energy levels are 
such that Ei < E 2 < E 3 . To couple the modes, it is sufficient to have two modulating 
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fields of the possible three that would be in resonance with the corresponding transition 
frequencies. In general, there can be three detunings A 21 , A 31 , and A 32 satisfying the 
resonance conditions (JBJ). For the transition amplitudes ()12)). we have j3i 2 , (3x3, and (3 23 . 
There are six amplitudes (fTTjl for with i ^ j. 

Realizing the coupling of three modes by two driving fields, we can create three types of 
resonant systems, which may be called, by analogy with the similar situations for resonant 
atoms, as cascade, V-type, and A-type systems jU^. In the cascade-type generation, the 
modes with the transition frequencies lu 2 i and a; 32 are coupled. In the V-type case, 
the transition frequencies are io 2 x and u)-&\. And for the A-type system, the transition 
frequencies are io%x and oj 32 . 

Substituting into the GPE ([!} the presentation (jl4|) and the corresponding driving 
fields (0), we employ the averaging technique [13]. This procedure is absolutely the same 
as has been thoroughly described in jUIE], so we do not repeat it here. The result is the 
system of equations for the coefficient functions 

l — =[OLl2\C2\ +«13|C 3 | C1+F1, 



dt 

UC 2 / I |2 I I o\ „ 

1 ~dt = I" 21 ' 01 ' +a23 l C3 l J C 2 + ^2, 

i = («3i|ci| 2 + a 32 |c 2 | 2 ) c 3 + F 3 , (36) 

where the terms Fi depend on the type of the generation method. For the cascade 
generation, we have 

F x = ~ (3 12 c 2 e 4A2lt , 
F2 = \ #2 ex + \ (3 23 c 3 e iA -* , 

Fz = \ P* 23 c 2 e- jA32 ' . (37) 

If we set here «i 3 , a 3 x, a 23 , a 32 , and /3 23 to zero, we return to the studied earlier two-mode 
case jHl E] • In the case of the V-type coupling, 

Fx = \ (3 12 c 2 e* A -' + \ (3 13 c 3 e lA -* , 
F 2 = \ P* 12 c x , F 3 = 1 /3* 3 cj e- jA -* . (38) 



And for the A-type generation, 

F x = \ Pis c 3 e iA ^ , F 2 = 1 /5 23 c 3 e iA ^ , 

F 3 = i cx + i /?* 3 c 2 e"* A ^ . (39) 

In addition, from the normalization of the function 1)140. we have 

|ci| 2 + |c 2 | 2 + |c 3 | 2 = l. (40) 
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Each Q = Cj(i) defines the dynamics of the corresponding fractional mode population 
lc-1 2 

It is important to stress that Eqs. (|36j) are obtained by employing the standard av- 
eraging method [HI]), taking into account the existence of two time scales, slow and fast, 
related to the inequalities (fT3"|) and (fTo'|). In this averaging procedure, one substitutes 
expansion ()14|) into the GPE (0) and averages over time fastly oscillating functions. As a 
result of this, the so-called time-reversed or counter-rotating terms vanish, so that in the 
right-hand side of Eqs. ()36j) there are no terms like c\c\ or c^cf. It is well known that, 
even if such terms would be added, they do not produce a significant change in the solu- 
tions. Omitting these terms is essentially equivalent to the widely known rotating-wave 
approximation [33]. The derivation of equations for c n (t), with all related mathematical 
details has been thoroughly described in Refs. P 15] l2~2"j. 

Though Eqs. (|36|) look differently for different types of mode generation related to 
distinct terms F^, the mathematical structure of these equations is, actually, the same. We 
may notice the following symmetry properties. The V^-type equations can be obtained 
from the cascade-type ones by interchanging the indices 1 and 2 and by replacing f3 2 i 
by (3{ 2 . Similarly, the A- type equations can be derived from the cascade- type equations 
by interchanging the indices 2 and 3, with the replacement f3 32 —> ^23 ■ The relation 
= — Aji has to be taken into account. Because of this symmetry, it is sufficient 
to consider just one type of Eqs. (jHSI), for instance, that corresponding to the cascade 
generation. 

The functions Cj(t) are complex- valued. Hence, Eqs. (J3T?|) present a system of six 
differential equations. However, it is possible to show that they define a four- dimensional 
dynamical system. To prove this, we involve the notation 

Cj = \cj\ exp(i7Tj) , (41) 

where iij = iTj(t) is a real- valued phase. Also, we write 

Pij = b,^ exp(i7 ii ) , b^ = |/%| . (42) 

Introduce the population differences 

s= |c 2 | 2 -| Cl | 2 , p = |c 3 | 2 — |c 2 | 2 (43) 

and the relative phases 

X = 7T 2 - 7Tl + 712 + A 2i t , V = 7T 3 - 7T 2 + 7 23 + A 32 t . (44) 

The fractional mode populations can be expressed through the variables as 

M 2 = \ (1 - 2s - p) , \c 2 \ 2 = - (1 + s - p) |c 3 | 2 = ^ (1 + s + 2p). 
Then Eqs. f!36|) for the cascade generation can be reduced to the system of four equations 

-77 = \ \A + s -p (b 23 \J 1 + s + 2p sin y - 2b 12 \Jl-2s-p sin x) , 
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dp 1 
dt ~ 3 



1 + s — p (bi2 \Jl — 2s — p sin x — 25 23 \J 1 + s + 2p sin y^j 



dx 36i2 s cos x 1 / 1 + s + 2p 

— = Otis + 0\p H . - - 6 23 W t— cos y + J 2 

a* 2a/(1 + s-p)(1-2s-p) 2 I ! + s-j) 



% r 36 2 3PCOS?/ 1 1 — 2s — p . s 

— = a 2 p + d 3 s^ = + 7T b i2W-r^ cosx + d) 4 , (45) 



in which 

1 . 1 . , 

«1 = - («12 + «13 + 2 "21 - «23) , «2 = - («32 + «31 + 2 «23 ~ "21) , 

Si = ^ («21 - «12 + 2 "13 - 2 «23) , *2 = A 2 l + - («12 - «21 + «13 - "23) , 

<^3 = ^ ("23 - "32 + 2 «31 - 2 «2l) , h = A 32 + - (« 23 ~ «32 + «21 ~ «3l) • 

Thus, the set of six equations ()36j) is really equivalent to a four-dimensional dynamical 

system 

We have numerically investigated the behaviour of solutions to Eqs. (jUJ) for various 
parameters ctj, 6y, and <5j. The latter parameter, playing the role of an effective detuning, 
was assumed to be small, 5i <C 1. Different initial conditions have been considered in 
the range — 1 < s < 1, — 1 < p < 1, < Xq < 2ir, < y < 2-n. For small <C «j, 
the solutions for the population differences s and p demonstrate a kind of nonlinear Rabi 
oscillations in the mode locked regime, when s(t) and p(t) do not cross the zero line, 
being always either above or below it, depending on initial conditions. This mode locked 
regime is the same as in the case of two modes, studied earlier (HI CH H3 |2U 122] • The 
difference with the two-mode case is that the Rabi-type oscillations look slightly more 
complicated, being quasiperiodic but not periodic. Increasing bij results in the increase 
of the oscillation amplitudes, and, after a critical value of b^, the mode unlocked regime 
arises, when either s(t), or p(t), or both of them, oscillate in the whole interval [-1,1]. 
To illustrate these two regimes, we present numerical calculations accomplished for equal 
parameters = a with the initial conditions s = —1, p = 0, x = y = 0. Figured] 
shows the mode locked regime, when s(t) < for all t > 0, and p(t) is also negative for 
almost all t. Figure |2] demonstrates the mode unlocked regime, when s(t) G [—1, 1], while 
pit) is yet almost always negative. Increasing further fry leads to the situation, when pit) 
starts also oscillating in the interval [-1,1]. However, at large b^, the temporal behaviour 
of solutions becomes quite unstable resembling chaotic motion. To better understand 
what happens, it is necessary to study the phase portrait of the dynamical system, that 
is, one has to find the fixed points and to analyze their stability. However, before turning 
to this task, we will compare the results of the averaging method described above with 
those of a direct numerical simulation of the GPE for several special cases. 

5 Direct numerical simulation of time evolution 



Although the averaging technique is a well grounded mathematical method j4oJ |36], it 
is interesting to show that the main effects, described above, also appear in a direct 



12 



numerical simulation of the time evolution described by Eq. ((TJ . One important reason 
is that the results of the averaging procedure are basically independent of the details of 
the trapping potential, while a direct numerical simulation is affected by it. We consider 
two cases: in an anharmonic potential the conditions for the validity of Eq. (}3l)|) can well 
be met and reasonable agreement can be expected. On the other hand, in a harmonic 
trap the averaging method may break down. We will demonstrate both situations using 
the numerical methods and parameter settings which are described in the Appendix. 

5.1 Anharmonic trap with linear driving field 

The case of an anharmonic trapping potential is most suited to demonstrate the benefits 
of the few-mode averaging method. We considered a BEC moving in a potential of the 
form Uqz a with Uq = 10" 32 J//xm 4 and a linear driving potential of the form V\ = (3\z (in 
the notation of Eq. (JHJ)) which was tuned to resonance with the <po fi transition so 
that uji = {Ei — E )/h ps 610/s. In absence of the driving potential the BEC can occupy 
several stationary nonlinear coherent modes, three of which are displayed in Fig. El The 
respective parameters appearing in Eq. © are found to be a 12 = 121.7s" 1 , a^ 3 = 46.3s" 1 , 
a 2 i = 144.2s" 1 , ct 23 = 88.9s" 1 , a 31 = 91.8s" 1 , and ct 32 = 111.9s" 1 . 

As predicted by the averaging method our simulation showed a critical behaviour: For 
f3\ below a certain value /3 C ~ 7.3 x 10" 33 J//xm the population transferred to the first 
excited state is bounded to be smaller than about 0.5. An example for this behaviour can 
be seen in Fig. 0Ja) which corresponds to a slightly subcritical value of (3i = 6.9 x 10" 33 
J//zm. For f3\ > P c this upper bound suddenly disappears and the population of the 
first excited states can take almost all values between and 1. This can be seen in 
Fig. E] a) which displays the time evolution for a slightly supercritical driving field with 
ft = 7.7 X 10" 33 J/fim. 

The corresponding predictions of the averaging technique can be seen in Figs |Ub) and 
E3b), respectively. Since the period and amplitude of the populations |cj| 2 are strongly 
varying around the critical value of f3i, we have chosen the values (3\ = 1.0 x 10" 32 J//im 
for 0] b) and (3\ = 1.02 x 10" 32 J//im for El b), which are close to the critical value as 
predicted by the averaging method. The found critical values of (3 C are in reasonable 
agreement, the difference being of about 25%, which is, actually, the accuracy of the 
averaging technique for the given set of parameters. 

5.2 Harmonic trap with anharmonic driving field 

This case serves as an example of how the breakdown of the conditions (|T3"j) and (|15|) 
results in wrong predictions of the mode expansion method. The case of a harmonic trap 
is a very special case with this respect: if, in absence of interaction, the transition between 
two neighbouring modes is resonantly driven, then this is also the case for a transition 
between any other neighbouring modes. It therefore is never possible to consider only a 
small number of states because other states are quickly populated as well. Only if one 
is sufficiently far away from resonance a few-mode model can be expected to work well, 
but then the transition rate between those modes is also low. These conclusions remain 
qualitatively also valid in the presence of interaction. 
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To demonstrate the breakdown of the mode expansion method we consider a harmonic 
potential of the form U z (z) = m u 2 z 2 /2 with oj z = 600/s, and a screened cubic driving 
field which, in the notation of Eq. (jUJ), is given by Vi(z) = faz 3 exp(—z 2 /w 2 ), i = 1, 2, with 
w = 4.9/im. We introduced the exponential screening since a purely cubic driving field 
plus a harmonic potential is unbounded from below. A linear driving field cannot be used 
because, according to the theorem given above, only leads to an oscillation of the BEC 
without changing its shape. The two driving frequencies are chosen to be at resonance 
with the ipo — > (fx transition, to\ = (E\ — E )/h « 536/s, and the (p% — > (p 2 transition, 
uj 2 = (E 2 -E x )/htt 568/s. 

The relative strength of the two driving fields was chosen to achieve (tpi,V\Lpo) = 
(<f2,V2 l Pi) by setting the appropriate values for (3\ and /3 2 . The results are shown in 
Fig. For weak driving fields, & = 0.7869 x 1(T 33 J//im 3 and p 2 = 0.344 x 10~ 33 
J//zm 3 , one achieves reasonable agreement with the predictions of the averaging procedure, 
see Fig. Et). For a strong driving force with f3i = 1.96734 x 10~ 33 J//mi 3 and f3 2 = 
0.859975 x 10 -33 J//im 3 the excitation of higher modes becomes significant, which is not 
surprising because the driving field provides more energy which can also excite higher 
levels. This behaviour can be seen in Fig. EJd) which also indicates that the time evolution 
deviates strongly from the predictions of the averaging method. 



6 Stationary Solutions and Stability Analysis 

To find out the stationary solutions for the dynamics of the three-mode case and to analyse 
the stability of these solutions, it is convenient to work with the variables 

£ = 1^1 (7 = 1,2,3) (46) 

and the relative phases (jUj). The fractional mode populations are expressed through the 
amplitudes (j4^)l as rij = /?. 

Using the variables fj, x, and y one can rewrite Eqs. ()36|) for the cascade generation 
in the form 

dfi 1 . 

-j- = - - 612/1 siM + - 6 2 3 /3 smy , 

dh 1, , . 
^ = -2 Wssm,, 



d X A „ f2 1 „ f2 1 /„ „. \ f2 1 1 fl „„„„ z, f'< 



A 2 1 - «21 fl + "12 fl + ("13 - "23) fi + 6l2 J t f r L COSX ~ b 23 777- COS V 



dt * x ZWi y ±0 '°' J6 " 2/J2 2/ 
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dy fi f 2 — f 2 

— = A 32 + («21 - "3l) fl - "32 fl + "23 fl + 6l2 7TT COSX + b 23 t , , 2 COS?/ . (47) 
« t ^/2 ^/2/3 

There are here yet too many parameters, because of which the analysis of Eqs. (|47jl is yet 
too complicated. To simplify the consideration, we may involve a realistic approximation, 
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when there is no detuning, the amplitudes due to the applied alternating fields, are 
taken to be the same, and the parameters <Xij are close to each other. That is, we set 

a = ay , b = , Aij = . (48) 

It is also convenient to measure time in units of a -1 . To return back to time units, one 
should replace t by at. Then Eqs. (J4T)) reduce to 

dfx b . 

—— = - to SIM, 

dt 2 J 

dh b b 

—— = fi sin x H — U sin y , 

dh b 

—— = t 2 smy , 

^ = / 2 -/ 1 +&^-cosx-6— cosy, 

d 4 = fl-f" + b^- cosx + b^-f cosy. (49) 
In addition, because of Eqs. (|3U|) and (JIBJ), there is the relation 

fx + fl + /! = !■ (50) 

This relation is automatically supported by Eqs. (jj^) as well as by Eqs. (|4*7jl . provided it 
is valid for the initial /j(0). Equations ()49|) are much easier to analyse than Eqs. Q47]). At 
the same time, the mathematical structure of these equations is similar, so the behaviour 
of their solutions should be close to each other. 

In Eqs. (I49j) . there is the sole dimensionless parameter b, defined in Eq. (148)) . Since 
the parameter a can be positive as well as negative, then b can, generally, be of both signs 
too. But Eqs. ()49|) enjoy a nice symmetry property, being invariant under the change 

b — > —b , x-^x±7r, ?/ — > y ± 7r . (51) 

Therefore, in what follows, it is sufficient to consider only the case of positive b > 0. 
There exists a special solution of Eqs. (14911 . for which 



fx = fs , x = -y (52) 

for all times t > 0. Then the problem reduces to a kind of a two- mode case described by 
the equations 

dfx b 

—— = - to SIM, 

dt 2 J 

dh , , . 
— = -b fx sin x , 
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This reduction, however, becomes possible as a result of the approximation when all 
interaction amplitudes cty = a are assumed to be equal. More precisely, the minimal 
requirements necessary for the existence of solution (J52j) are 

«12 = «32 , «13 = «31 , &12 = &23 , A 2 1 = A 32 . (53) 

It is feasible, in principle, to choose such modulating fields and a trapping potential that 
Eqs. (JH3|) be valid. In addition, the initial conditions are to be such that /i(0) = /3(C)). 

The stationary solutions of Eqs. (|49jl are obtained by equating to zero their right- 
hand sides, keeping in mind that h 7^ and / 2 is not identically zero. The corresponding 
fixed-point equations for the phase differences are 

sinrc = siny = , cosx = ±1 , cosy = ±1 . (54) 

The equations for the amplitudes are rather cumbersome, and we shall not write them 
down. We shall solve these equations numerically, calculating the fixed-point amplitudes 
f* = f*(b) as functions of the pumping parameter b and then finding the related stationary 
solutions 

<{b) = |/*(&)| 2 (55) 

for the fractional mode populations. Simultaneously, we accomplish the stability analysis 
by calculating the Jacobian matrix for Eqs. (|4*9~j) . evaluated at the related fixed points. 
The eigenvalues of this matrix define the characteristic exponents characterizing the type 
of stability |47j . 

The first stationary solution is given by x* = y* = and n* = n%. The related mode 
populations ()55j) are shown in Fig. [7| This fixed point is neutrally stable, being a center 
with the characteristic exponents Ai = and imaginary A 2 = A3, A4 = A£. When b 
varies in the interval 0.1 < b < 1, the absolute values of A 2j 3 and A 4j 5 change in the range 
0.192 < |A 2j3 | < 0.865 and 0.351 < |A 4 , 5 | < 1-755. 

The second fixed point is defined by x* = y* = n and the branch I in Fig. |Sk) for 
n\ = rig and in Fig. (Hb) for n 2 . This fixed point is also a center with the characteristic 
exponents A x = and imaginary A 2 = A3, A 4 = AJj. For b G [0,1], one has |A 23 | ~ 
|A 4 , 5 | ~ 1. 

The third fixed point corresponds to x* = y* = n and the branch II in Figs. IHt) and 
IHb) for n\ = and n 2 , respectively. The characteristic exponents are Ai = and real 
A 2 = —A 3 , A 4 = — A 5 , which shows that this point is unstable. The absolute values are 
|A 2 , 3 | ~ |A 4)5 | ~ 0.1. The solution exists for < b < b* , with b* = 0.198431. 

The fourth fixed point is given by x* = y* = tc and the branch III in Figs. |Ht) and 
IHb) for n\ = n% and n 2 - The characteristic exponents are Ai = 0, imaginary A 2 = A3, 
and real A 4 = — A5, with |A 2j 3| ~ |A 4) s| ~ 0.1. This solution is unstable and exists under 
< b < b* , with 6* = 0.198431. 

The fifth fixed point is described by x* = y* = it and the branch I in Fig. for nl(b), 
n 2 (b), and n^(b), respectively. The characteristic exponents are A 4 = and imaginary 
A 2 = A3, A 4 = A5, which shows that this point is a center. The solution exists for b in 
the range < b < b* c , with b* = 0.639448, where |A 2)3 | ~ 1 and |A 4)5 ] < 0.967. 

The sixth fixed point corresponds to x* = y* = tt and the branch II in Fig. Elfor nl(b), 
n 2 (fc), and nl(b). The related characteristic exponents are A 4 = 0, imaginary A 2 = A3, and 
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real A4 = — A5. This tells that this point is unstable. The solution exists for < b < b*, 
with b* c = 0.639448, where |A 2)3 | < 0.896 and |A 4j5 | < 0.457. 

The stationary solutions for the case n* 7^ are, actually, invariant under the in- 
terchange of the indices 1 and 3. For concreteness, we accept that n\ > n\, as is shown 
in Fig. Formally, there is also the seventh fixed point, for which x* = y* = 2tt and 
n\ = n\. But this is just the first point shifted in x* and y* by 2ti. 

Selecting from all available fixed points only those that corresponds to stable stationary 
solutions, we have three such cases: the point x* = y* = 0, n\ = 72,3, depicted in Fig. [7] the 
point x* = y* = it, n\ = rig, shown in Fig. [HI as the branch I; and the point x* = y* = tt, 
n\ > rig, presented as the branch I in Fig. El 

Recall that these three stable fixed points exist under the validity of conditions ([53)1 . In 
reality, it can be quite difficult to satisfy these conditions exactly. But if these conditions 
are not valid, then the sole stable point that remains is the point x* = y* = n, n* > n^, 
which is shown as the branch I in Fig. This stationary solution exists only for b < b*. 
For larger pumping parameters b > b*, there are no stable (or neutrally stable) fixed 
points. Hence, the motion for such large b will be chaotic. 

7 Harmonic Generation and Parametric Conversion 

Employing the presentation (fT4*j) in the GPE and involving the averaging technique 
|4~3] . we come to the equations for c n (t), like Eqs. ()3b|). This procedure defines an initial 
approximation for c n (t), which can be called the guiding centers and labelled as cffl(t). 
It is possible to obtain corrections to the guiding centers by using the following steps 
of the averaging technique [13 ES] • Then we can find the higher approximations (t) 
describing the fractional mode populations |c^(t)| 2 . To obtain the higher approximations 
for c n (t), we may proceed as follows. 

Let us present the solution to the temporal equation (JTJ) as 



where again c n (t) is a slow function of time, compared to the fast function ip n (r,t). For 
any fast function of time F(t), we define the averaging 



(56) 
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(57) 



over fast oscillations. In particular, 




Ul ^ UJ 2 ■ 



The functions ip n (r,t) can be taken such that 




(58) 



Then the amplitudes c n (t) in Eq. can be obtained as 




(59) 
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and the normalization 

J2K(t)\ 2 = i (60) 

n 

is valid. For instance, taking tp n (r,t) in the form 

^\r,t) = <p n (r) exp (- % - E n tj , (61) 

we come back to Eq. (|14j) . with all Eqs. (j58|l to (j60|) being evidently satisfied. 

Equation can be employed as a relation for an iterative procedure defined by the 
rule 

c { n k+1) (t) = <J ^>(r,t) <p<*>(r,f) dv > , (62) 

where 

^(r,t)^c^(t)^(v,t). (63) 

n 

Starting from the guiding centers c$>(t) and the form (JBTj) . we get 

c£ D W = cf)(t), (64) 

which follows from the rule (|f)2jl . 

To derive the second-order approximation for c n (t), we write 

?P(r,t) = K(r) + X n(r,t)] exp (- 1 E n . (65) 

Here <p n (r) is a stationary topological mode given by the eigenproblem (jlj), and x n (r, t) 
has to be found by substituting Eq. (JoT)|) into the GPE ((TJ). In the case of modulating 
field ©, the correcting term x n (r,t) can be written as 



Xn(r. 



^EM'l^ + W^] > ( 66 ) 



where j = 1,2, . . ., and all c^- > can be ordered so that < u\ < u>2 < • . .. The 
functions u n j(r) and f n j(r) satisfy the equations 



[H[f n ] - E n + N A s \p n \ 2 - tkUjJ u nj + N A s ip 2 n v nj = - - cp n B* , 

(6[<p n ] -E n + N A s \p n \ 2 + htoj) v nj + N A s { Vn f u n j = --(PnBj, (67) 

where <p n = <p n (r). Using the first-order form in Eqs. (J53|) and (JHH), we obtain the 
second-order approximation 

c?(t) = c^(t) + £ cg)(t) [(^ Mrm ) + «, ^ m )] Afa - u nm ) + 

m { i 

H Kv 3 ^ v* mi ) + (u ni , y? m )] A(^ + u nm ) + 

i 
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+ H [( U ni, U m j) A(^j - UJj + UJ nm ) + (u ni , V mj )A(uJi + UJj + UJ nm ) + 
ij 

+ C*4> U mj )A(Ui + UJj - UJ nm ) + (?4, V* mj ) A((Ji ~ UJj - UJnm)] } , (68) 

in which (u, v) denotes the corresponding scalar product. This formula is valid for an 
arbitrary number of alternating resonant fields. In the case of just one resonant field, there 
should be no summation over the indices % and j in Eq. ([68)1 . Formula (|68J) shows what 
are the conditions on the alternating fields, which induce resonant transitions between 
topological modes. These conditions depend on the number of the modulating fields 
involved. 

Consider the case of two resonantly coupled topological modes, when 

cf(t) = (n^l,2), (69) 

and is nonzero only for n — 1, 2, as follows from the averaging technique for the guiding 
centers [OJIEI- And suppose that there is the sole alternating field with a frequency uj\ = to. 
Then the nontrivial population amplitudes are 

cP(t) = c f\t) + cf\t) [(ip u v*) + K,¥>2)] A(o; - uj 21 ) + 4 0) (t)( Ml ,^)A(2cu - uj 21 ) , 

#(*) = 4 0) (t) + cS 0) (t) [(<p 2 ,ui) + (yln)] A(uj-uj 21 ) + cf\t)(v;, Ul )A(2uj-uj 21 ) , (70) 

which results from Eq. 1)68)1 . taking into account that CJ12 = — ct>2i- 

If there are two alternating fields, so that j = 1,2 in Eq. ()68j) . then for the nonzero 
population amplitudes, one has 

cP = c f\t) + 4 0) {[(^ 2 \) + (u luV2 )] A( Wl - cu 21 )+ 

+ [(^1^22) + ( M i2 ? ^2)] A(u; 2 - W21) + (un, v 21 )A(2u 1 - a; 2 i) + (ui2, v 22 )A(2uj 2 - ^21) + 
+ [(u n , v 22 ) + (ui2, v 21 )] A(ui +uj 2 - uj 21 ) + [(un, U21) + (v 22 , v u )] A(uj 2 -u x - uj 2 i)} , 

= c { 2 \t) + cf ] {[(ip 2 ,u n ) + A(wi - cu 2 i) + 

+ [(^2, W12) + (^22) <Pi)] A(u; 2 - w 2 i) + («2i> «n)A(2wi - cj 2 i) + «i2)A(2u;2 - w 2 i)+ 

+ [(«2i, W12) + (^22> A(^l + UJ 2 - U21) + [{U21, U12) + Ki, «2a)] A(^ 2 - Wi - w 2 i)} • 

(71) 

Expressions 1)70)1 show that one alternating field, with a frequency u;, can induce 
transitions between two topological modes, provided that one of the following equations 
is valid: 

uj = uj 2 i , 2uj = uj 2 i . (72) 

And Eqs. (|71)1 demonstrate that two alternating fields, with frequencies uj\ and uj 2 , can 
realize intermode transitions under one of the resonance conditions: 

lO\ = UJ 2 \ , UJ 2 = UJ 2 \ , 

2u)\ = uj 21 , 2uj 2 = uj 2 i , 
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UJi + UJ 2 = U>21 > Co> 2 — Co"i = U21 ■ (73) 

Going to the higher-order approximations in the iterative procedure (j62j). we see that, in 
addition to the standard resonance conditions as lj = u 2 \ or uj,i = u 2 i, there appear the 
conditions of harmonic generation 

ku = u 2l (k = 2,3,...) , (74) 

when there is the sole alternating field, or the conditions of parametric conversion 

J2( ±UJ j) = ^21 , (75) 

3 

if several alternating fields modulate the trapping potential. 

The prediction of nonlinear harmonic generation using the averaging method is in 
excellent agreement with a respective numerical simulation for a harmonic potential of 
the form U z (z) = m uj 2 z 2 /2, with u z = 600/s. To demonstrate harmonic generation we 
employed a harmonic driving field of the form (jUJ) with V\(z) = (3z 2 and driving frequency 
u>i = (E 2 —Ei)/(2h) m 552/s. For a sufficiently strong driving field ((3 = l.lxlCr 32 J//xm 2 , 
roughly half the strength of the trapping potential) we observed a strongly enhanced 
population of the second excited mode (excitation of the first excited state is forbidden 
by parity conservation). The result for the time averaged coefficients is shown in Fig. ITU1 
and demonstrates a population of up to 30 % in this mode, far more than is to be expected 
in absence of interaction. As in the case of an anharmonic driving field (see above) we 
also find a substantial (30 %) excitation of other modes for a strong driving field. 

The terms harmonic generation and parametric conversion are used here by analogy 
with nonlinear optics, where there exist analogous phenomena [13]. Similar effects occur 
as well for elementary excitations of Bose-Einstein condensates |481 14*§1 15D] . As we have 
shown, such effects may also arise for topological coherent modes. 

In conclusion, the resonant generation of topological modes seems to provide a feasible 
mechanism for creating novel states of trapped Bose gases, containing the combinations 
of several such modes. Here, we have considered the resonant generation realized by 
modulating the trapping potential. Note that another possibility of imposing resonant 
perturbations could be done by inducing periodic variations in the atomic scattering 
length jHIl E2] • The latter case requires that the atoms are in a state close to a Feshbach 
resonance. 

As we mentioned in the Introduction, a dipole topological mode was successfully gen- 
erated in experiment [23] . The resonant method of vortex creation was analysed in detail 
in Thus, the possibility of generating a single topological mode is well justified. 

The mode can be of any nature, whether dipole or vortex. The method of resonant gener- 
ation works for any mode. Our numerical simulations for the GPE, with realistic physical 
parameters, show that the simultaneous generation of several topological modes is also 
feasible. 
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A Appendix. Description of numerical algorithms 

Our numerical simulations are performed for a BEC that obeys the GPE with an external 
potential of the form U (r) = Uj_(r±)+U z (z), i.e., a transverse and a longitudinal trapping 
potential. We consider a quasi one-dimensional situation for which the energy required to 
create excited states of the transverse trapping potential U±(t±) is much higher than any 
other energy. The wavefunction then takes the approximate form (p(r, t) = <p±(r±) ip(z, t), 
where <p±(r±) is the transverse ground state. By integrating over the transverse coordi- 
nates r_|_ one can show that ip fulfills a one- dimensional GPE of the form 



with an effective coupling parameter A z = A s J |<^±(r_i_)| 4 dr±, whereby the integral over 
the transverse ground state is of the order of the transverse width squared. The details 
of this reduction can be found in review [2|. 

In this appendix we shall present the methods used for a direct numerical solution of 
the one-dimensional GPE for different situations discussed in the preceding sections. We 
generally consider a BEC of 100 87 Rb atoms (mo = 1.45 x 10 -25 kg, a s = 5.4 nm) with a 
transverse width of 7 /iin. Our ID simulation described a time evolution of 500 ms using 
10 7 discrete time steps. It was performed on a grid of 512 spatial points extending over a 
range of 26 /xm. 

A.l Numerical simulation of time evolution and determination 
of nonlinear coherent modes 

To simulate the time evolution governed by Eq. (J7BJ) we have employed the well-known 
split-step Fourier method [SB]- As initial condition we used the ground-state nonlinear 
coherent mode (po(z) associated with the potential U z (z) in absence of a driving field 
V(z,t). This ground state as well as higher nonlinear coherent modes were numerically 
determined by an imaginary-time propagation followed by a self-consistent field (SCF) 
method. 

The imaginary-time propagation was described, for instance, in Ref . j^l] • Very roughly 
speaking it is a method to obtain a kind of a small-temperature (corresponding to large 
imaginary time) collective wavefunction which approaches the ground state wavefunction 
for very large times. It usually gives excellent results for the ground state but fails to 
reproduce excited nonlinear coherent modes unless certain symmetry requirements are 
fulfilled. We used this method for intermediate imaginary times to obtain a better trial 
wavefunction ip(°) for the following SCF procedure. This step is not necessary but helps 
to improve the speed of the numerical algorithm. 

The SCF method consists in inserting the trial wavefunction ip^ into the nonlinear 
interaction part of Eq. (|76|). It thus provides a kind of potential term and turns Eq. 
(J76j) into a linear equation for ip. We then solved for the stationary eigenstates of this 
equation. One of these states is then chosen to become a new trial wavefunction cp^ 
and then the procedure is iterated until the overlap between cp^ and y?( n+1 ) deviates by 
less than e from 1 (we have chosen e = 10~ 14 ). The choice of the appropriate new trial 




(76) 
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wavefunction depends on the problem at hand. In the case of a ID Schrodinger equation 
with a confining potential, the eigenstates are real and the number N Q of zeros of each 
nonlinear coherent mode turned out to be a good criterion for the next trial wavefunction. 
To determine the ground state we picked the N = eigenstate in each iteration. For the 
excited coherent modes N = 1 and 2 was the respective choice. 

The convergence of the SCF method is the slower the higher the excited state one is 
looking for, and it even can fail to converge for No > 2. We therefore implemented a 
convergence acceleration scheme based on Anderson mixing as described in Ref. In 
short, this scheme enhances the usual SCF algorithm by adding a kind of memory in the 
sense that y?( n+1 ) may depend not only on (p( n > but on previous trial wavefunctions, too. In 
practice, this memory lasts for 3-6 iterations. To implement this scheme one first writes 
(p{n+i) as a superposition of the previous trial wavefunctions. Then one minimizes the 
deviation of this function from the previous trial wavefunction with respect to the super- 
position coefficients. This leads to a simple linear system of equations for the coefficients 
of the superposition. Anderson mixing can speed up convergence by orders of magnitude 
and turned out to be very suitable for our case. The three nonlinear coherent modes 
of lowest energy for a harmonic potential and an anharmonic potential are displayed in 
Fig. El 

We have checked if the nonlinear coherent modes obtained in this way are indeed 
stationary solutions of the nonlinear Schrodinger equation by propagating them for a 
sufficiently long time using the split step algorithm. All modes were found to be perfectly 
stationary: their density does not change and their phase remains spatially homogeneous 
apart from the sharp phase jump by ir that appears at the zeros of the excited states, and 
apart from areas with extremely low densities (less than about 10~ 13 of the peak density) 
where numerical errors lead to a strongly fluctuating phase. 

A. 2 Comparison to averaging procedure 

To compare the result of the direct numerical simulation with that obtained by the aver- 
aging procedure we first projected the numerical solution ip(z, t) found using the split-step 
method on the respective nonlinear coherent mode to obtain as an intermediate result the 
non-averaged coefficients Cj(tj) := ((pi,ip(ti)) exp(iEiti/h) with ti = lAt,l = 0,1,2,.... 
Since the nonlinear coherent modes (fi are not orthogonal the sum Y^i \ci{ti)\ 2 is not unity 
and the individual coefficients are oscillating on short time scales. 

To obtain the coefficients Cj(tj) corresponding to the averaging procedure we stored 
the wavefunction at a total number of N t time steps (N t was typically on the order of 
1000) and defined the averaged coefficients as 

I N a 

Ci(T n ):=— $>(T n + f,), (77) 
iVa i=i 

where N a is the number of steps we averaged over and T n := N a nAt, n = 0,1,2,.... 
Generally N a = 10 to N a = 32 resulted in a fairly smooth time evolution of the coefficients 
Ci(t) which can directly be compared to the results of the averaging procedure. 
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Figure 1: Mode locked regime for the case of three coupled nonlinear modes. Parameters 
are = a, b,^ = 0.35a, and 5i = 0. Initial conditions are s = —1, Po = 0, x = y = 0. 
Time is measured in units of a -1 . Shown are the population differences s(t) (solid line) 
and p(t) (dashed line). 




Figure 2: Mode unlocked regime for the three-mode case. All parameters and initial 
conditions are the same as in Fig. 1, except by = 0.55a. Again, time is shown in units of 
or 1 . The population differences: (a) s(t); (b) p(t). 
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Figure 3: The three nonlinear coherent modes of lowest energy for 100 87 Rb atoms 
in two different trapping potentials. The solid line corresponds to the ground state, the 
dashed (dotted) line to the first (second) excited mode, respectively, a) harmonic trap 
with frequency uo z = 600/s, b) anharmonic trap of the form UqZ 4 with Uq = 10~ 32 J//im 4 . 




Figure 4: Mode coefficients averaged over a time of 3 ms for a BEC in an anharmonic 
trap oc z 4 driven resonantly by a slightly subcritical linear potential, a) results of direct 
numerical simulation, b) predictions of the averaging method (see text for more details). 
Solid line: |c (t)| 2 , dashed line: |ci(t)| 2 , dotted line: |c 2 (t)| 2 . 
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Figure 6: Mode coefficients averaged over a time of 5.3 ms for a BEC in a harmonic 
trap driven resonantly by a screened cubic potential, a) weak driving force with fa = 
0.7869 x KT 33 J//im 3 and (3 2 = 0.344 x 10~ 33 J//im 3 . b) strong driving force with 
A = 1.96734 x 10- 33 J//im 3 and f3 2 = 0.859975 x lO" 33 J/^m 3 . Solid line: |c (t)| 2 , dashed 
line: |ci(t)| 2 , dotted line: |c 2 (t)| 2 . 




Figure 7: Stationary solutions for the case x* = y* = and n\ = as functions of the 
pumping parameter b: stable n\(b) = nl(b) (solid line), stable nl(b) (dashed line). 
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a) 



b) 



Figure 8: Stationary solutions related to the fixed point x* = y* = it and n\ = n% as 
functions of the pumping parameter b. Stable branch I (solid line), unstable branch II 
(dashed line), and unstable branch III (dashed line): (a) n\{b) = n%{b); (b) ri^ib). 




Figure 9: Stationary solutions associated with the fixed point x* = y* = it and n\ > n| 
as functions of b. Stable branch I (solid line) and unstable branch II (dashed line): (a) 
nl(b); (b) n$(6); (c) n$(6). 




Figure 10: Time evolution of a BEC trapped and driven by harmonic potentials, a) mode 
coefficients averaged over a time of 5.3 ms (the lines have the same meaning as in Fig. 0J), 
b) spatial probability density at different time steps. The effect of harmonic generation 
of population in the second excited state can clearly be seen in both figures. 
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